##
## Jun Koga Sudduth 
##
## JPR article "Coup risk, coup-proofing and leader survival"
## jun.koga@strath.ac.uk
##
## Models in Tables I and II in the main text 

library(R2WinBUGS)
library(MCMCpack)
library(foreign)
library(MASS)
library(coda)
library(AER) # ivreg
library(plm)
library(zoo)
library(ggplot)
 


##
## Coup Risk Model (Table I in the main text)
##
##

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus")

dat <- read.table(file="Data/dat12132015_outsample.txt")

names(dat)


# Set up to estimate in BUGS;
dat <- dat[dat$year>=1965 & dat$year<=2003,]

dim(dat)  # 5927 observations


# missing data
dat <- subset(dat,  is.na( lag_grgdpch_gled)==F  & is.na(lag1regime)==F & is.na(anycoupyrs)==F)

attach(dat)
sum(is.na(lag1regime))



# N
N <- dim(dat)[1]  # 5927 observations
N

#Coup Attempt Variable

Anycoup <- dat$any_coup  # any coup


# Lag 1 Log of GDP per capita

lag1econGled <- dat$lln_rgdpch_gle
  
# democracy

lagDemocracy <- dat$lag1democracy

lagMilitary  <- dat$lag1mildic

# years after coups

Coupyrs <- dat$anycoupyrs
 

##
## For Prediction
##

# Key Economy variable
hist(lag1econGled, breaks=30, col="blue")
summary(lag1econGled)
xz1 <- seq(5, 11, 0.01)
length(xz1)# 601

mean(GrowthGled, na.rm=TRUE) #  0.0190378
mean(Coupyrs, na.rm=TRUE)  # 17.56414


# Key variables (economy) range and mean/median of other controls

xz1 <- cbind( 0.019,  xz1, 0, 1, 17)






##
## Couprisk 1
##

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/R")

inits <- function() { list(d1=rnorm(1),
                           d2=rnorm(1),
                           d3=rnorm(1),
                           d4=rnorm(1),
                           d5=rnorm(1)
                           )}


model.file <- "JPR_R2_couprisk_base.txt"

bugs.model.couprisk1 <- bugs(data=c("Anycoup", "lag1econGled", "lagDemocracy", "lagMilitary", "Coupyrs", "N"),
                    inits=inits,
                    parameters.to.save=c("d", "q"),    # "q", "y.hat"
                    model.file=model.file,
                    n.chains=3,
                    n.iter=3000,
                    bugs.directory="C:/Program Files (x86)/WinBUGS14/",
                    debug=TRUE)

plot(bugs.model.couprisk1)

print(bugs.model.couprisk1)


# Post-WinBugs Analysis

dat$qmean1 <- bugs.model.couprisk1$mean$q # Note that you need to first save 'estimated coup risk" in the Data. Then you'll make "q".
write.dta(dat,file="C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data/dat_R2_base.dta")

##
## Results
##

# 95%
meand <- bugs.model.couprisk1$mean$d
sdd <- bugs.model.couprisk1$sd$d

lowd <- NULL
highd <- NULL
for (i in 1:5){
l <- quantile(bugs.model.couprisk1$sims.list$d[,i],0.025)
h <- quantile(bugs.model.couprisk1$sims.list$d[,i], 0.975)
lowd <- c(lowd, l)
highd <- c(highd, h)
}

data.frame(meand, sdd , lowd, highd)
#                             meand         sdd      lowd     highd
#                              0.83286585 0.661304585 -0.429580  2.164500
# lag GDP per capita           -0.40400471 0.089861034 -0.584450 -0.229450
# Democracy                      0.07128797 0.182812300 -0.275210  0.430140
# Military                         0.38200372 0.151381765  0.083247  0.680030
# Coup Yrs                       -0.07838684 0.009113505 -0.097348 -0.061768

# 90%
meand <- bugs.model.couprisk1$mean$d
sdd <- bugs.model.couprisk1$sd$d

lowd <- NULL
highd <- NULL
for (i in 1:5){
l <- quantile(bugs.model.couprisk1$sims.list$d[,i],0.05)
h <- quantile(bugs.model.couprisk1$sims.list$d[,i], 0.95)
lowd <- c(lowd, l)
highd <- c(highd, h)
}

#         meand         sdd      lowd     highd
#1  0.83286585 0.661304585 -0.296160  2.027400
#2 -0.40400471 0.089861034 -0.568140 -0.251620
#3  0.07128797 0.182812300 -0.217680  0.390480
#4  0.38200372 0.151381765  0.116060  0.621960
#5 -0.07838684 0.009113505 -0.094102 -0.064332







##
## Paramilitary DV (Model 1 in Table II in the main text)
##

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data")

dat <- read.dta(file="C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data/dat_R2_base.dta")

# Set up to estimate in BUGS;

dat <- dat[dat$year>=1968 & dat$year<=2003,] # Restrict the year where paramilitary data is available (1968-2003)
dim(dat)  # 5927 observations

## Missing Data

dat <- subset(dat,  is.na( lag1qmean1  )==F & is.na(lag_priowar  )==F & is.na( lag_priocivil )==F & is.na(  lag1pararatio )==F & is.na(lag1regime)==F   )


N <- dim(dat)[1]  # 5927 observations

# DV
qlogPara <- qlogis(dat$pararatio  )   # qlogis(p)=logit(p)=log(p/1-p)=log(p)-log(1-p)  where p is paramilitary ratio here.
lag1qlogPara <-  qlogis( dat$lag1pararatio )

# democracy

lagDemocracy <- dat$lag1democracy

# war
lagPrioCivil <- dat$lag_priocivil
lagPrioWar <- dat$lag_priowar

# Coup risk
lagq1 <- dat$lag1qmean1



## For Prediction
##

xp <- seq(0,0.25, .001) # coup risk range for a prediction
length(xp)   # 251

summary(lag1qlogPara)

xp <- cbind(xp, 0, 0, 0,  -0.62020)# No war/ no democracy

 
#  Model
#

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/R")


inits <- function() { list(
                           b1=rnorm(1),
                           b2=rnorm(1),
                           b3=rnorm(1),
                           b4=rnorm(1),
                           b5=rnorm(1) ,
                           b6=rnorm(1),
                           tau=1
                           )}

model.file <- "JPR-12212015-para1_prediction.txt"

bugs.model.para1 <- bugs(data=c("qlogPara", "lagq1", "lagPrioWar", "lagPrioCivil", "lagDemocracy", "lag1qlogPara", "xp", "N"),
                    inits=inits,
                    parameters.to.save=c("b", "mu", "y.hat"),
                    model.file=model.file,
                    n.chains=3,
                    n.iter=1000, # note you need more than 10000 to get enough convergence.
                    bugs.directory="C:/Program Files (x86)/WinBUGS14/",
                    debug=TRUE)


plot(bugs.model.para1)


print(bugs.model.para1)


# 95% ok for Lag1 Coup risk

meanb  <- bugs.model.para1$mean$b
sdb <- bugs.model.para1$sd$b

lowb <- NULL
highb <- NULL
for (i in 1:6){
l <- quantile(bugs.model.para1$sims.list$b[,i],0.025)
h <- quantile(bugs.model.para1$sims.list$b[,i], 0.975)
lowb <- c(lowb, l)
highb<- c(highb, h)
}

data.frame(meanb, sdb , lowb, highb)
#                            meanb         sdb        lowb       highb
#                        -0.02065743 0.020935736 -0.06133575  0.02131725
# Lag Coup Risk          -0.81586063 0.309662167 -1.44057500 -0.23004500
# Lag Int War            -0.05035739 0.057604462 -0.16191000  0.06169200
# Lag Civil War           0.02596166 0.030785900 -0.03408200  0.08386700
# Democracy              -0.09145107 0.026586496 -0.14465250 -0.04141350
# Lag DV                 0.84691053 0.009961788  0.82754250  0.86645250
#

 
    


# 95% ok for Lag1 Coup risk

meanb  <- bugs.model.para1$mean$b
sdb <- bugs.model.para1$sd$b

lowb <- NULL
highb <- NULL
for (i in 1:6){
l <- quantile(bugs.model.para1$sims.list$b[,i],0.05)
h <- quantile(bugs.model.para1$sims.list$b[,i], 0.95)
lowb <- c(lowb, l)
highb<- c(highb, h)
}

data.frame(meanb, sdb , lowb, highb)

#    meanb         sdb       lowb      highb
#1 -0.02065743 0.020935736 -0.0555435  0.0131360
#2 -0.81586063 0.309662167 -1.3380500 -0.3209250
#3 -0.05035739 0.057604462 -0.1444200  0.0424805
#4  0.02596166 0.030785900 -0.0248725  0.0759245
#5 -0.09145107 0.026586496 -0.1364150 -0.0479895
#6  0.84691053 0.009961788  0.8313000  0.8631000
>







##
## Counterbalancing DV (Model 5 in Table II in the main text)
##

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data")

 dat <- read.dta(file="C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data/dat_R2_base.dta")
# Set up to estimate in BUGS;

dat <- dat[dat$year>=1968 & dat$year<=2003,] # Restrict the year where paramilitary data is available (1968-2003)
dim(dat)  # 5927 observations

## Missing Data

dat <- subset(dat, is.na( lag1qmean1   )==F & is.na(lag_priowar  )==F & is.na( lag_priocivil )==F & is.na( lag1cp_BS  )==F & is.na(lag1regime)==F   )


N <- dim(dat)[1]  # 5927 observations


# DV
CBalance <- dat$cp_BS
lag1CBalance <- dat$lag1cp_BS

# democracy

lagDemocracy <- dat$lag1democracy

# war
lagPrioCivil <- dat$lag_priocivil
lagPrioWar <- dat$lag_priowar

# Coup risk
lagq1 <- dat$lag1qmean1
lag2q1 <- dat$lag2qmean1


## For Prediction
##

xp <- seq(0,0.25, .001) # coup risk range for a prediction
length(xp)   # 251

summary(lag1CBalance)

xp <- cbind(xp, 0, 0, 0, 0.01081 )# No war/ no democracy


# Model

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/R")


inits <- function() { list(
                           b1=rnorm(1),
                           b2=rnorm(1),
                           b3=rnorm(1),
                           b4=rnorm(1),
                           b5=rnorm(1) ,
                           b6=rnorm(1),
                           tau=1
                           )}

model.file <- "JPR-12212015-counter1_prediction.txt"

bugs.model.balance1 <- bugs(data=c("CBalance", "lagq1", "lagPrioWar", "lagPrioCivil", "lagDemocracy", "lag1CBalance", "xp", "N"),
                    inits=inits,
                    parameters.to.save=c("b", "mu" , "y.hat"),
                    model.file=model.file,
                    n.chains=3,
                    n.iter=1000,
                    bugs.directory="C:/Program Files (x86)/WinBUGS14/",
                    debug=TRUE)


plot(bugs.model.balance1)
print(bugs.model.balance1)





# 95% ok for Lag1 Coup risk

meanb  <- bugs.model.balance1$mean$b
sdb <- bugs.model.balance1$sd$b

lowb <- NULL
highb <- NULL
for (i in 1:6){
l <- quantile(bugs.model.balance1$sims.list$b[,i],0.025)
h <- quantile(bugs.model.balance1$sims.list$b[,i], 0.975)
lowb <- c(lowb, l)
highb<- c(highb, h)
}

data.frame(meanb, sdb , lowb, highb)

> data.frame(meanb, sdb , lowb, highb)
#                                  meanb         sdb        lowb       highb
#                              0.06827849 0.024037798  0.02180950  0.11691500
# Lag Coup Risk                -0.74844410 0.354496883 -1.45957500 -0.07547575
# Lag Int War                   0.04826219 0.066828352 -0.08117675  0.17826250
# Lag Civil War                 0.09332600 0.035440359  0.02417775  0.15985250
# Lag Democracy                 -0.09659496 0.030525512 -0.15765750 -0.03871850
# Lag DV                         0.89348400 0.008946094  0.87589500  0.91125250
>

 



# 90% ok for Lag1 Coup risk

meanb  <- bugs.model.balance1$mean$b
sdb <- bugs.model.balance1$sd$b

lowb <- NULL
highb <- NULL
for (i in 1:6){
l <- quantile(bugs.model.balance1$sims.list$b[,i],0.05)
h <- quantile(bugs.model.balance1$sims.list$b[,i], 0.95)
lowb <- c(lowb, l)
highb<- c(highb, h)
}

data.frame(meanb, sdb , lowb, highb)
#      meanb         sdb      lowb     highb
#1  0.06827849 0.024037798  0.028170  0.108500
#2 -0.74844410 0.354496883 -1.358000 -0.176055
#3  0.04826219 0.066828352 -0.060873  0.155915
#4  0.09332600 0.035440359  0.034946  0.150805
#5 -0.09659496 0.030525512 -0.148600 -0.046376
#6  0.89348400 0.008946094  0.879300  0.908500
>








##
##
## Prediction Plot
##
##

# The Effect of Coup Risk on Paramilitary

windows(15,8)
par(mfrow=c(1,2), mar=c(4,4,1.2,0.1))

xp <- seq(0, 0.25, .001)
length(xp) # 251


plot(xp, plogis(bugs.model.para1$mean$y.hat) , xlim=c(0,0.26),ylim=c(0.3, 0.4), main="Paramilitary", xlab="Coup Risk at t-1", ylab="Predicted Paramilitary Proportion at t", col="white")
lines(xp, plogis(bugs.model.para1$mean$y.hat),   col="blue")

ly.hat<- NULL
hy.hat <- NULL
for (i in 1:251){
l <- quantile(bugs.model.para1$sims.list$y.hat[,i],0.025)
h <- quantile(bugs.model.para1$sims.list$y.hat[,i],0.975)
ly.hat <- c(ly.hat, l)
hy.hat <- c(hy.hat, h)
}
lines(xp, plogis(ly.hat), lwd=1, lty=3, col="blue")
lines(xp, plogis(hy.hat), lwd=1, lty=3, col="blue")

legend(0.15, 0.4, c("Paramilitary Proportion", "2.5% CI", "97.5% CI"), lwd=1, lty=c(1,3,3), col="blue", cex=.85)



##
## The Effect of Coup Risk on Counterbalancing

xp <- seq(0, 0.25, .001)
length(xp) # 251

#windows()
#par(mar=c(5,4,1,1))

plot(xp, bugs.model.balance1$mean$y.hat , xlim=c(0,0.26), ylim=c(-0.2, 0.2),main="Counterbalancing", xlab="Coup Risk at t-1", ylab="Predicted Counterbalancing at t", col="white")
lines(xp, bugs.model.balance1$mean$y.hat,   col="blue")

ly.hat<- NULL
hy.hat <- NULL
for (i in 1:251){
l <- quantile(bugs.model.balance1$sims.list$y.hat[,i],0.025)
h <- quantile(bugs.model.balance1$sims.list$y.hat[,i],0.975)
ly.hat <- c(ly.hat, l)
hy.hat <- c(hy.hat, h)
}
lines(xp, ly.hat, lwd=1, lty=3, col="blue")
lines(xp, hy.hat, lwd=1, lty=3, col="blue")

legend(0.15, 0.18, c("Counterbalancing", "2.5% CI", "97.5% CI"), lwd=1, lty=c(1,3,3), col="blue", cex=.85)



# Save
savePlot(filename="C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Paper/main_main", type = "pdf")







 

##
## Model 2 in Table II in the main text (nonlinear model) 
##

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data")
 
dat <- read.dta(file="C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data/dat_R2_base.dta")


# Set up to estimate in BUGS;

dat <- dat[dat$year>=1968 & dat$year<=2003,] # Restrict the year where paramilitary data is available (1968-2003)
dim(dat)  # 5927 observations

## Missing Data

dat <- subset(dat, is.na( lag1qmean1  )==F & is.na(lag_priowar )==F & is.na(lag_priocivil)==F     &   is.na(  lag1pararatio )==F & is.na(lag1regime)==F   )


N <- dim(dat)[1]  # 5927 observations

# DV
qlogPara <- qlogis(dat$pararatio  )   # qlogis(p)=logit(p)=log(p/1-p)=log(p)-log(1-p)  where p is paramilitary ratio here.
lag1qlogPara <-  qlogis( dat$lag1pararatio )

# democracy

lagDemocracy <- dat$lag1democracy
 
# war
lagPrioCivil <- dat$lag_priocivil
lagPrioWar <- dat$lag_priowar

# Coup risk
lagq1 <- dat$lag1qmean1

Sqrlagq1 <- (dat$lag1qmean1)^2

#
# Paramilitary Model
#

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/R")


inits <- function() { list(
                           b1=rnorm(1),
                           b2=rnorm(1),
                           b3=rnorm(1),
                           b4=rnorm(1),
                           b5=rnorm(1) ,
                           b6=rnorm(1),
                           b7=rnorm(1),
                           tau=1
                           )}

model.file <- "JPR-nonlinear-para1a.txt"

bugs.model.para1 <- bugs(data=c("qlogPara", "lagq1", "Sqrlagq1", "lagPrioWar", "lagPrioCivil", "lagDemocracy", "lag1qlogPara", "N"),
                    inits=inits,
                    parameters.to.save=c("b", "mu"),
                    model.file=model.file,
                    n.chains=3,
                    n.iter=1000,  
                    bugs.directory="C:/Program Files (x86)/WinBUGS14/",
                    debug=TRUE)


plot(bugs.model.para1)

print(bugs.model.para1)

 

# 95% ok for Lag1 Coup risk

meanb  <- bugs.model.para1$mean$b
sdb <- bugs.model.para1$sd$b


lowb <- NULL
highb <- NULL
for (i in 1:7){
l <- quantile(bugs.model.para1$sims.list$b[,i],0.025)
h <- quantile(bugs.model.para1$sims.list$b[,i], 0.975)
lowb <- c(lowb, l)
highb<- c(highb, h)
}

data.frame(meanb, sdb , lowb, highb)

#                              meanb         sdb        lowb       highb
#                              -0.02001190 0.023634886  -0.0645855  0.02775025
# lag Coup Risk                  -0.77247557 0.757637228  -2.2040500  0.65171500
# Coup Risk ^2                   -0.41454337 4.937710568 -10.4210000  9.40890000
# lag Int War                    -0.05185683 0.059937715  -0.1707725  0.06850300
# lag Civil War                    0.02729453 0.030990330  -0.0315065  0.08640850
# lag Democracy                   -0.09213728 0.026196871  -0.1417150 -0.04139425
# lag DV                         0.84749140 0.009734856   0.8286900  0.86660000
>


# 90% ok for Lag1 Coup risk

meanb  <- bugs.model.para1$mean$b
sdb <- bugs.model.para1$sd$b


lowb <- NULL
highb <- NULL
for (i in 1:7){
l <- quantile(bugs.model.para1$sims.list$b[,i],0.05)
h <- quantile(bugs.model.para1$sims.list$b[,i], 0.75)
lowb <- c(lowb, l)
highb<- c(highb, h)
}

data.frame(meanb, sdb , lowb, highb)
#               meanb         sdb       lowb     highb
#1 -0.02001190 0.023634886 -0.0585615 -0.004210
#2 -0.77247557 0.757637228 -2.0290500 -0.229625
#3 -0.41454337 4.937710568 -8.6828000  2.799500
#4 -0.05185683 0.059937715 -0.1493150 -0.011495
#5  0.02729453 0.030990330 -0.0219820  0.048065
#6 -0.09213728 0.026196871 -0.1361000 -0.073720
#7  0.84749140 0.009734856  0.8314000  0.853700





##
## Model 3 in Table II in the main text (ParaDV controlling lag 1 Coup attempt)
##

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data")

 dat <- read.dta(file="C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data/dat_R2_base.dta")

#
##
## Results   (from 1965- 2
# Set up to estimate in BUGS;

dat <- dat[dat$year>=1968 & dat$year<=2003,] # Restrict the year where paramilitary data is available (1968-2003)
dim(dat)  # 5927 observations

## Missing Data

dat <- subset(dat, is.na(lag1any_coup ) ==F & is.na( lag1qmean1  )==F & is.na(lag_priowar  )==F & is.na( lag_priocivil )==F & is.na(  lag1pararatio )==F & is.na(lag1regime)==F   )


N <- dim(dat)[1]  # 5927 observations

# DV
qlogPara <- qlogis(dat$pararatio  )   # qlogis(p)=logit(p)=log(p/1-p)=log(p)-log(1-p)  where p is paramilitary ratio here.
lag1qlogPara <-  qlogis( dat$lag1pararatio )

# lag Coup Attempt
lagCoupAttempt <- dat$lag1any_coup


# democracy

lagDemocracy <- dat$lag1democracy

 

# war
lagPrioCivil <- dat$lag_priocivil
lagPrioWar <- dat$lag_priowar

# Coup risk
lagq1 <- dat$lag1qmean1

#
# Paramilitary Model
#

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/R")


inits <- function() { list(
                           b1=rnorm(1),
                           b2=rnorm(1),
                           b3=rnorm(1),
                           b4=rnorm(1),
                           b5=rnorm(1) ,
                           b6=rnorm(1),
b7=rnorm(1),
                           tau=1
                           )}

model.file <- "JPR-12212015-para1_attempt.txt"

bugs.model.para1 <- bugs(data=c("lagCoupAttempt",    "qlogPara", "lagq1", "lagPrioWar", "lagPrioCivil", "lagDemocracy", "lag1qlogPara", "N"),
                    inits=inits,
                    parameters.to.save=c("b"),
                    model.file=model.file,
                    n.chains=3,
                    n.iter=1000,  
                    bugs.directory="C:/Program Files (x86)/WinBUGS14/",
                    debug=TRUE)


plot(bugs.model.para1)


print(bugs.model.para1)



# 95% ok for Lag1 Coup risk

meanb  <- bugs.model.para1$mean$b
sdb <- bugs.model.para1$sd$b

lowb <- NULL
highb <- NULL
for (i in 1:7){
l <- quantile(bugs.model.para1$sims.list$b[,i],0.025)
h <- quantile(bugs.model.para1$sims.list$b[,i], 0.975)
lowb <- c(lowb, l)
highb<- c(highb, h)
}

data.frame(meanb, sdb , lowb, highb)

#                                 meanb         sdb        lowb       highb
#                              -0.01947868 0.021008668 -0.06045250  0.02063000
# lag1  Coup Risk               -0.86377661 0.306134634 -1.45452500 -0.25944250
# lag1 Int War                   -0.05053440 0.059893637 -0.16937250  0.06973775
# lag1 Civil War                  0.02300854 0.031446721 -0.03651700  0.08346525
# lag1 Democracy                 -0.09081059 0.026548689 -0.14420000 -0.03861400
# lag1 DV                        0.84719073 0.009707974  0.82764750  0.86615750
# lag1 Coup Attempt                 0.04822569 0.066728687 -0.07676175  0.17645250
>



# 90% ok for Lag1 Coup risk

meanb  <- bugs.model.para1$mean$b
sdb <- bugs.model.para1$sd$b

lowb <- NULL
highb <- NULL
for (i in 1:7){
l <- quantile(bugs.model.para1$sims.list$b[,i],0.05)
h <- quantile(bugs.model.para1$sims.list$b[,i], 0.95)
lowb <- c(lowb, l)
highb<- c(highb, h)
}

data.frame(meanb, sdb , lowb, highb)
    meanb         sdb       lowb      highb
1 -0.01947868 0.021008668 -0.0545105  0.0141420
2 -0.86377661 0.306134634 -1.3660500 -0.3538100
3 -0.05053440 0.059893637 -0.1479150  0.0457395
4  0.02300854 0.031446721 -0.0294375  0.0762815
5 -0.09081059 0.026548689 -0.1335000 -0.0470600
6  0.84719073 0.009707974  0.8310000  0.8628000
7  0.04822569 0.066728687 -0.0592500  0.1563050
>












##
## Model 4 in Table II in the main text (Para DV , Excluding lag1 Coup Attempt ==1)
##

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data")

 dat <- read.dta(file="C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data/dat_R2_base.dta")

#
# Set up to estimate in BUGS;

dat <- dat[dat$year>=1968 & dat$year<=2003,] # Restrict the year where paramilitary data is available (1968-2003)
dim(dat)  # 5927 observations

## Missing Data

dat <- subset(dat, is.na( lag1qmean1  )==F & is.na(lag_priowar  )==F & is.na( lag_priocivil )==F & is.na(  lag1pararatio )==F & is.na(lag1regime)==F   )
dat <- subset(dat, lag1any_coup==0  )

N <- dim(dat)[1]  # 5927 observations

# DV
qlogPara <- qlogis(dat$pararatio  )   # qlogis(p)=logit(p)=log(p/1-p)=log(p)-log(1-p)  where p is paramilitary ratio here.
lag1qlogPara <-  qlogis( dat$lag1pararatio )

# democracy

lagDemocracy <- dat$lag1democracy



# war
lagPrioCivil <- dat$lag_priocivil
lagPrioWar <- dat$lag_priowar

# Coup risk
lagq1 <- dat$lag1qmean1

#
# Paramilitary Model
#

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/R")


inits <- function() { list(
                           b1=rnorm(1),
                           b2=rnorm(1),
                           b3=rnorm(1),
                           b4=rnorm(1),
                           b5=rnorm(1) ,
                           b6=rnorm(1),
                           tau=1
                           )}

model.file <- "JPR-12212015-para1.txt"

bugs.model.para1 <- bugs(data=c("qlogPara", "lagq1", "lagPrioWar", "lagPrioCivil", "lagDemocracy", "lag1qlogPara", "N"),
                    inits=inits,
                    parameters.to.save=c("b"),
                    model.file=model.file,
                    n.chains=3,
                    n.iter=1000,
                    bugs.directory="C:/Program Files (x86)/WinBUGS14/",
                    debug=TRUE)


plot(bugs.model.para1)


print(bugs.model.para1)


# 95% ok for Lag1 Coup risk

meanb  <- bugs.model.para1$mean$b
sdb <- bugs.model.para1$sd$b

lowb <- NULL
highb <- NULL
for (i in 1:6){
l <- quantile(bugs.model.para1$sims.list$b[,i],0.025)
h <- quantile(bugs.model.para1$sims.list$b[,i], 0.975)
lowb <- c(lowb, l)
highb<- c(highb, h)
}

data.frame(meanb, sdb , lowb, highb)

#                            meanb        sdb        lowb       highb
#                            -0.01516967 0.02159223 -0.05665675  0.02607525
#  Coup Risk                  -0.87325794 0.32077449 -1.50305000 -0.25894250
# lag Int War                -0.05614598 0.05678397 -0.17085250  0.05548375
# lag Civil War               0.01108599 0.03276983 -0.05283200  0.07657925
# lag Democracy               -0.09362449 0.02636163 -0.14430000 -0.04107975
# lag DV                       0.84834927 0.00961061  0.83014750  0.86665250



# 90% ok for Lag1 Coup risk

meanb  <- bugs.model.para1$mean$b
sdb <- bugs.model.para1$sd$b

lowb <- NULL
highb <- NULL
for (i in 1:6){
l <- quantile(bugs.model.para1$sims.list$b[,i],0.05)
h <- quantile(bugs.model.para1$sims.list$b[,i], 0.95)
lowb <- c(lowb, l)
highb<- c(highb, h)
}

data.frame(meanb, sdb , lowb, highb)

#             meanb        sdb       lowb      highb
#1 -0.01516967 0.02159223 -0.0506690  0.0197095
#2 -0.87325794 0.32077449 -1.4100000 -0.3600650
#3 -0.05614598 0.05678397 -0.1491150  0.0346160
#4  0.01108599 0.03276983 -0.0436905  0.0662470
#5 -0.09362449 0.02636163 -0.1374150 -0.0502620
#6  0.84834927 0.00961061  0.8331900  0.8637100
>











##
## Model 6 in Table II in the main text (Counterbalancing DV nonlinear)
##

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data")

 dat <- read.dta(file="C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data/dat_R2_base.dta")


# Set up to estimate in BUGS;

dat <- dat[dat$year>=1968 & dat$year<=2003,] # Restrict the year where paramilitary data is available (1968-2003)
dim(dat)  # 5927 observations

## Missing Data

dat <- subset(dat, is.na( lag1qmean1  )==F & is.na(lag_priocivil )==F & is.na(lag_priowar )==F   &   is.na( lag1cp_BS)==F & is.na(lag1regime)==F   )

N <- dim(dat)[1]  # 5927 observations


# DV
CBalance <- dat$cp_BS
lag1CBalance <- dat$lag1cp_BS

# democracy

lagDemocracy <- dat$lag1democracy

# war
lagPrioCivil <- dat$lag_priocivil
lagPrioWar <- dat$lag_priowar

# Coup risk
lagq1 <- dat$lag1qmean1

Sqrlagq1 <- (dat$lag1qmean1)^2

#
#

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/R")


inits <- function() { list(
                           b1=rnorm(1),
                           b2=rnorm(1),
                           b3=rnorm(1),
                           b4=rnorm(1),
                           b5=rnorm(1) ,
                           b6=rnorm(1),
                           tau=1
                           )}

model.file <- "JPR-nonlinear-counter1a.txt"

bugs.model.counter1 <- bugs(data=c("CBalance", "lagq1",   "Sqrlagq1", "lagPrioWar", "lagPrioCivil",  "lagDemocracy", "lag1CBalance", "N"),
                    inits=inits,
                    parameters.to.save=c("b", "mu"),
                    model.file=model.file,
                    n.chains=3,
                    n.iter=1000,  
                    bugs.directory="C:/Program Files (x86)/WinBUGS14/",
                    debug=TRUE)


plot(bugs.model.counter1)
print( bugs.model.counter1)
 

# 95% ok for Lag2 Coup risk

meanb  <- bugs.model.counter1$mean$b
sdb <- bugs.model.counter1$sd$b

lowb <- NULL
highb <- NULL
for (i in 1:7){
l <- quantile(bugs.model.counter1$sims.list$b[,i],0.025)
h <- quantile(bugs.model.counter1$sims.list$b[,i], 0.975)
lowb <- c(lowb, l)
highb<- c(highb, h)
}

data.frame(meanb, sdb , lowb, highb)

>
#                                       meanb        sdb         lowb      highb
#                                    0.07104780 0.02709050   0.02035325  0.1267625
#Lag Coup Risk                       -0.85935692 0.84658832  -2.46705000  0.7320550
#Coup Risk ^2                          0.67560372 5.47943748 -10.36000000 11.5905000
#Civil War                             0.09453157 0.03565831   0.02728725  0.1625100
#Int War                                0.04723856 0.06955629  -0.09080075  0.1869050
#Democracy                              -0.09737050 0.03009840  -0.15435250 -0.0394025
# lag DV                                  0.89390867 0.00872767   0.87720000  0.9105525
>



# 90% ok for Lag2 Coup risk

meanb  <- bugs.model.counter1$mean$b
sdb <- bugs.model.counter1$sd$b

lowb <- NULL
highb <- NULL
for (i in 1:7){
l <- quantile(bugs.model.counter1$sims.list$b[,i],0.05)
h <- quantile(bugs.model.counter1$sims.list$b[,i], 0.95)
lowb <- c(lowb, l)
highb<- c(highb, h)
}

data.frame(meanb, sdb , lowb, highb)

#   meanb        sdb       lowb     highb
#1  0.07104780 0.02709050  0.0275685  0.116705
#2 -0.85935692 0.84658832 -2.2631000  0.547325
#3  0.67560372 5.47943748 -8.6136500  9.218200
#4  0.09453157 0.03565831  0.0378560  0.152310
#5  0.04723856 0.06955629 -0.0658195  0.159015
#6 -0.09737050 0.03009840 -0.1476200 -0.047438
#7  0.89390867 0.00872767  0.8796000  0.908005




##
## Model 7 in Table II in the main text 
## (Counterbalancing DV controling lag1 coup attempt )
##

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data")
 
 dat <- read.dta(file="C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data/dat_R2_base.dta")

# Set up to estimate in BUGS;

dat <- dat[dat$year>=1968 & dat$year<=2003,] # Restrict the year where paramilitary data is available (1968-2003)
dim(dat)  # 5927 observations

## Missing Data

dat <- subset(dat, is.na( lag1any_coup  )==F &  is.na( lag1qmean1   )==F & is.na(lag_priowar  )==F & is.na( lag_priocivil )==F & is.na( lag1cp_BS  )==F & is.na(lag1regime)==F   )
#dat <- subset(dat, lag1any_coup==0 & any_coup==0)

N <- dim(dat)[1]  # 5927 observations


# DV
CBalance <- dat$cp_BS
lag1CBalance <- dat$lag1cp_BS

# democracy

lagDemocracy <- dat$lag1democracy

# war
lagPrioCivil <- dat$lag_priocivil
lagPrioWar <- dat$lag_priowar

# Coup risk
lagq1 <- dat$lag1qmean1
lag2q1 <- dat$lag2qmean1

lagCoupAttempt <- dat$lag1any_coup


#
#

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/R")


inits <- function() { list(
                           b1=rnorm(1),
                           b2=rnorm(1),
                           b3=rnorm(1),
                           b4=rnorm(1),
                           b5=rnorm(1) ,
                           b6=rnorm(1),
                           b7=rnorm(1),
                           tau=1
                           )}

model.file <- "JPR-12212015-counter1_attempt.txt"

bugs.model.para1 <- bugs(data=c("lagCoupAttempt",   "CBalance", "lagq1", "lagPrioWar", "lagPrioCivil", "lagDemocracy", "lag1CBalance", "N"),
                    inits=inits,
                    parameters.to.save=c("b" ),
                    model.file=model.file,
                    n.chains=3,
                    n.iter=1000,  
                    bugs.directory="C:/Program Files (x86)/WinBUGS14/",
                    debug=TRUE)


plot(bugs.model.para1)
print(bugs.model.para1)

  

# 95% ok for Lag1 Coup risk

meanb  <- bugs.model.para1$mean$b
sdb <- bugs.model.para1$sd$b

lowb <- NULL
highb <- NULL
for (i in 1:7){
l <- quantile(bugs.model.para1$sims.list$b[,i],0.025)
h <- quantile(bugs.model.para1$sims.list$b[,i], 0.975)
lowb <- c(lowb, l)
highb<- c(highb, h)
}

data.frame(meanb, sdb , lowb, highb)

#
#                                  meanb         sdb        lowb       highb
#                                 0.06945755 0.024008587  0.02266400  0.11680000
# lag1 Coup Risk                  -0.80494204 0.350311840 -1.48352500 -0.12413500
# lag1 Int War                     0.04804785 0.069482491 -0.08984075  0.18755250
# lag1 Civil War                   0.08981379 0.036190300  0.02168250  0.15945750
# lag1 Democracy                   -0.09584722 0.030465759 -0.15705250 -0.03551175
# lag1 DV                         0.89372440 0.008758855  0.87649500  0.91125250
 # lag1 Coup Attempt                  0.05672120 0.076758367 -0.08702400  0.20425750
>




# 90% ok for Lag1 Coup risk

meanb  <- bugs.model.para1$mean$b
sdb <- bugs.model.para1$sd$b

lowb <- NULL
highb <- NULL
for (i in 1:7){
l <- quantile(bugs.model.para1$sims.list$b[,i],0.05)
h <- quantile(bugs.model.para1$sims.list$b[,i], 0.95)
lowb <- c(lowb, l)
highb<- c(highb, h)
}

data.frame(meanb, sdb , lowb, highb)

#  meanb         sdb       lowb     highb
#1  0.06945755 0.024008587  0.0289985  0.107705
#2 -0.80494204 0.350311840 -1.3790500 -0.223640
#3  0.04804785 0.069482491 -0.0648805  0.159710
#4  0.08981379 0.036190300  0.0297135  0.150925
#5 -0.09584722 0.030465759 -0.1444050 -0.045814
#6  0.89372440 0.008758855  0.8794000  0.908000
#7  0.05672120 0.076758367 -0.0668610  0.181010
>











##
## Model 8 in Table II in the main text (Counterbalancing DV excluding lag1 coup attempt ==1)
##

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data")

 dat <- read.dta(file="C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/Data/dat_R2_base.dta")


# Set up to estimate in BUGS;

dat <- dat[dat$year>=1968 & dat$year<=2003,] # Restrict the year where paramilitary data is available (1968-2003)
dim(dat)  # 5927 observations

## Missing Data

dat <- subset(dat, is.na( lag1qmean1   )==F & is.na(lag_priowar  )==F & is.na( lag_priocivil )==F & is.na( lag1cp_BS  )==F & is.na(lag1regime)==F   )
dat <- subset(dat, lag1any_coup==0)

N <- dim(dat)[1]  # 5927 observations


# DV
CBalance <- dat$cp_BS
lag1CBalance <- dat$lag1cp_BS

# democracy

lagDemocracy <- dat$lag1democracy

# war
lagPrioCivil <- dat$lag_priocivil
lagPrioWar <- dat$lag_priowar

# Coup risk
lagq1 <- dat$lag1qmean1
lag2q1 <- dat$lag2qmean1

#
#

setwd("C:/Users/Jun/Dropbox/Dissertation/Bayesian_prospectus/R")


inits <- function() { list(
                           b1=rnorm(1),
                           b2=rnorm(1),
                           b3=rnorm(1),
                           b4=rnorm(1),
                           b5=rnorm(1) ,
                           b6=rnorm(1),
                           tau=1
                           )}

model.file <- "JPR-12212015-counter1.txt"

bugs.model.para1 <- bugs(data=c("CBalance", "lagq1", "lagPrioWar", "lagPrioCivil", "lagDemocracy", "lag1CBalance", "N"),
                    inits=inits,
                    parameters.to.save=c("b" ),
                    model.file=model.file,
                    n.chains=3,
                    n.iter=1000,  
                    bugs.directory="C:/Program Files (x86)/WinBUGS14/",
                    debug=TRUE)


plot(bugs.model.para1)
print(bugs.model.para1)

 
# 95% ok for Lag1 Coup risk

meanb  <- bugs.model.para1$mean$b
sdb <- bugs.model.para1$sd$b

lowb <- NULL
highb <- NULL
for (i in 1:6){
l <- quantile(bugs.model.para1$sims.list$b[,i],0.025)
h <- quantile(bugs.model.para1$sims.list$b[,i], 0.975)
lowb <- c(lowb, l)
highb<- c(highb, h)
}

data.frame(meanb, sdb , lowb, highb)
#                          meanb         sdb        lowb       highb
#                            0.07292348 0.024787428  0.02731375  0.12257250
# lag Coup Risk             -0.85675329 0.366576629 -1.58557500 -0.15661250
# Lag Int War              0.03963651 0.065754988 -0.09318625  0.16887250
# Lag Civil War               0.07893546 0.037661947  0.00519955  0.15415750
# lag Demoracy                -0.09478361 0.030156291 -0.15300000 -0.03434475
# lag DV                       0.89631373 0.008576293  0.88000000  0.91235750
>



# 90% ok for Lag1 Coup risk

meanb  <- bugs.model.para1$mean$b
sdb <- bugs.model.para1$sd$b

lowb <- NULL
highb <- NULL
for (i in 1:6){
l <- quantile(bugs.model.para1$sims.list$b[,i],0.05)
h <- quantile(bugs.model.para1$sims.list$b[,i], 0.95)
lowb <- c(lowb, l)
highb<- c(highb, h)
}

data.frame(meanb, sdb , lowb, highb)

#  meanb         sdb       lowb      highb
#1  0.07292348 0.024787428  0.0324425  0.1127100
#2 -0.85675329 0.366576629 -1.4690000 -0.2730850
#3  0.03963651 0.065754988 -0.0680725  0.1447200
#4  0.07893546 0.037661947  0.0164575  0.1421100
#5 -0.09478361 0.030156291 -0.1445200 -0.0450865
#6  0.89631373 0.008576293  0.8823000  0.9101100





